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Abstract 

We map Eigen model of biological evolution [Naturwissenschaften 58, 465 (1971)] into a one- 
dimensional quantum spin model with non-Hermitean Hamiltonian. Based on such a connection, 
we derive exact relaxation periods for the Eigen model to approach static energy landscape from 
various initial conditions. We also study a simple case of dynamic fitness function. 
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Eigen model of asexual evolution is one of the main mathematical models in this 

field. In this model individuals have offsprings, that are subjected to mutation that connects 
with a selection rule. In his original work Eigen found an error threshold similar to the 
critical point in critical phenomena such that when the mutation is larger than the error 
threshold the organism can not survive. Later, statistical mechanics has been applied to 
investigate the discrete time version of the original model . Franz and Peliti I5| derived 
another important result in the Eigen model: concentration of individuals around the peak 
configuration. 

In the parallel mutation-selection model, an alternative to the Eigen model, a mutation 
mechanism and a selection mechanisms are two independent processes that take place con- 
currently (J. Baake et al. []| proved that for the parallel mutation-selection scheme, the time 
evolution equation for the frequencies of different species is equivalent to the Schrodinger 
equation in imaginary time for quantum spins in a transverse magnetic field. Based on 
such a connection, recently we used Suzuki- Trotter formalism jq] to study both statics and 
dynamics of the model with a single peak fitness function In the present Letter, we will 
extend such study to the Eigen model JlJ by reexpressing the Eigen model's dynamics via 
quantum chain problem, then solving the dynamics to obtain exact relaxation periods for 



the Eigen model. 



environments 
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he dynamic aspects play important role during the evolution in changing 
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12|. Thus such aspects in the Eigen model have been considered in 



recent works Il4j , in which approximate formulas for the relaxation periods have been 
found and applied to describe a virus-immune system coevolution. Our equations for exact 
relaxation periods are consistent with approximate formulas in Refs. ^| for the case of 
one mutation per replication. 

As in Ref. jlj, the genome configuration is specified by a sequence of N spin values 
Sfc = ±1, 1 < k < N. We denote the i-th genome configuration by Si = (si, s 2 , sjy) and 
the probability of the i-th genome at time t is given by = Pi(t) and the fitness is 
the average number of offspring's per unit time. In our language, the chosen fitness is a 
function / that operates on the genome configuration Si, i.e., = f(Si). 

In the Eigen model, elements of the mutation matrix represent the probability that 
an offspring produced by state j changes to state i, and the evolution is given by the set of 
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2 N coupled equations for 2 N probabilities pi 

dp- 2N 2N 

3=1 3=1 

Here p t satisfies E?=iP* = 1 and Q tj = q N - d (^)(l - q ) d ^) w i t h d(i,j) = (N - Ejli 4 s j)/ 2 
being the Hamming distance between Si and Sj. The parameter 1 — q describes the efficiency 
of mutations. For the parallel mutation-selection model, the dynamics is given by 

dp- 2N 2N 

3=1 3=1 

where rriij are the elements of the mutation matrix = 70 for d(i,j) = 1, TOy = —Njq for 
i = j, and = for d(i,j) > 1. 

Eigen found that it is enough to solve Eq. (JJJ for only linear parts Let us decompose 
the first, linear part of Eq. $1} via mutations to the fixed length d(i,j) = I: 

dp- N 

-^ = J2 J2 Q&iPi- ( 3 ) 
1=0 j,d(i,j)=i 

The second sum is over all configurations having Hamming distance I from the peak config- 
uration. Using the relation }2 i=1 Qij = 1, we can show that when pi satisfies Eq. J3J), then 



ms W) (4) 

satisfies Eq. (JIJ. We can compare Eq. (j3J) with Eq. (J2J) without the last nonlinear term. 
The terms I = 1 and I = in Eq. (jSJ) correspond, respectively, to the first and second terms 
in Eq. (J2J). In Eq. (jSJ), there are terms with higher level I > 2 spin flips. Baake et al. []| 
mapped Eq. (J2J) into a model of quantum spin chain. Here we will use the same method to 
map the model of Eqs. (0) and © into a quantum spin model with additional higher level 
spin flip terms. 

Let us reformulate the system of Eq. Q. As we identify configuration Sj with a collection 
of spins s^-.s^r — il and define fitness function / as Tj — /(s|..s^) = /(Sj). Let us consider 
vectors in the Hilbert space of iV quantum Pauli spins. With the pi of Eq.(jHJ), we connect a 
vector in Hilbert space zli=iPi\^i > - Then r 3 - — > f(al..a z N ). The / spin flip term in Eq. 
(JHJ) can be identified with a matrix element < Sj\Di\Si > of quantum operator 

Dl = q N - d ^\l-q) d ^ Yl <-<■ ( 5 ) 

l<ii<..tj<JV 
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Thus Eq. (3) is equivalent to Scrodinger equation: 

N _ 2^ 



N 2" 



Q 

1=1 (l<h<i 2 ..it<N) 

2 N 2 N 

j tZ ZvAm>=-H z ZvAm> (6) 



and Eq.(4) to: 

Z = £ < S l \e~ m \S j >p\ 



o 



Y. ] <S l \e- Hi \S 3 >p" 3 
Pi = — ^ , (7) 

where a denotes the spin operator and \S > is the standard notation for the spin state. One 
can multiply Eq. (6) from the left by < Si\ and obtain Eq. (J3J). 
For the single-peaked fitness function, we take 

f(S 1 ) = A, and /($) = 1 for 1, (8) 

with Si = (+1,-1-1, ...,+1), which is equivalent to choosing 



f ( Sl ) = i+(A-i)fi^r (9) 

at the limit p — > oo. A careful look at the Hamiltonian of Eq. (jBJ) reveals that it is non- 
Hermitean. But we will mainly work with the matrix elements between Si ^ S\ and Sj ^ Si 
and for these situations we can miss the multiplier f(af..a^f) = 1. For that sector of Hilbert 
space, Hamiltonian is Hermitean. To investigate the dynamics, we are using the matrix 
elements of Hamiltonian 

- < Si\H\Si >= Ae" 7 ; 
< Si\H\Sj >=< Si\Hdiff\Sj >, i 1; 

N 2 N 

-% / = ^ + ^e-^) f £ , (10) 

where / is identity operator, 7 = — TV ln(g) « iV(l — q) N. For us only terms Z <C iV are 
relevant, therefore the substitution ^[(1 — q)/q] 1 — > e^^/N) 1 can be applied. 

To calculate matrix elements of T(t) = e _m , one should introduce the Suzuki- Trotter 



formalism 



To perform analytical calculation, it is more convenient to use Eq. Q 



d ^Xi(t)\Si >= -H diff J2%i(tM > (12) 

i=l i=2 



for the fitness function and Eq. ()10jl . For any value of p an exact method of Suzuki- 
Trotter formalism ^| can map the system to the problem in classical statistical mechanics. 
Moreover, for the large values of p it is well known that the problem is drastically simplified. 
For the quantum p-spin interactions in a transverse magnetic field, Goldschmidt |15J] has 
found that all the order parameters (magnetizations) are either 1 or and one should take 
either only transverse interaction terms {af^.af ) or only the longitudinal one ( e~ 7 [l + (A — 
a i/N) p ])- Therefore, we can work with system of Eq. (jlOj) using the following trick. 
With exponential accuracy of order 1/2 N , it is possible to neglect the erf terms in Eq. © 
and get 

< Sxle'^Sx >~ exp[(Ar 7 )t]. (11) 

Matrix elements < Si\e~ Ht \Sj > for i ^ 1 can be replaced with exponential accuracy by 
< Si\exp[—Hdifft\\Sj >. Equation 

dt 

is equivalent to Eq. (J3j) with rj = 1 for j = 2, . . . , 2 and r\ = 0. Then we derive that 

2 jv 2^ 

^Xi(t) = exp[t] ^Xi. (13) 

i=2 i=2 

From Eqs. (fTT| and Eq. (|l3jl . we have p\ ~ exp[(Ae _7 )t] and ^2 i=2 Pi ~ e*. Therefore, we 
derive the Eigen's exact formulae for the error threshold, 

A > e 7 . (14) 

Let us calculate now the transition probabilities < Sj\exp[—Hdifft]\Si > between two 
states with the total number of M flips between configurations Si = {s\..s l N } and Sj = 
{sl-.s^} and define m = 1 — 2M/N. We will show below that the model can be solved at 

1 ~ (1 - q) « 1. (15) 

For the finite (1 — m), we guess that the relaxation time t is of order N and define 

T = ie- 7 /iV. (16) 

There are N(l + m)/2 spins without flips (+1 spins) and N(l — m)/2 flipped spins (-1 
spins). Let us denote by hi the term of I spin flip in the Hamiltonian. To calculate the 



matrix element < Sj\exp[—Hdifft]\Si >=< Sj\ exp[— t ^ hi]]\Si >, let us use an equality 
exp[aof of^ . . . of] = cosh[a][l + tanhfajer^cr? . . . of] and expand the product keeping terms 
till the M-th degree: 

- M\ 

< S j \e- ,H '"'\S, >« J2 £ J^T 

K=l h+..l K =M 

cosh( 7 T)-tanh( 7 T)^ fl ^ } (17) 

We find via the saddle point the principal term in the expression of Eq. (17) among all 
distributions with different Zj. We keep cosh,tanh only for the one spin flip terms. We 
calculate also the combinatorics of insertion into M site box combination of l\ single points, 
Z 2 duplets, ...Ik k plets, which satisfy the constraint 

M 

J2ih = M. (18) 

i=l 

We can take the constraint of Eq. (fTH|) into account via a Lagrange parameter A and write 
li as XiN . For the logarithm of a typical term for summation in Eq. (17), we have 

N(f)(T, m, 7) = AT [In cosh(7T) + x\ ln(tanh(7T)) 

1 — m 1 — m 1— m /rr ,s N 

+ ^— ln — ^ g ln(x^!/T) - s<) 

8=2 

A/ 1 - 

+ ln7^ zxj - Si lnxi + x\ + A(^ke, ; — )]. (19) 



i=2 

The extremum conditions for Xi of Eq. ()19j) give: 



xi = tanh( 7 T)z/7, ib* = TV, i > 2, (20) 



where z = r ye x . Using formulas: J2i^2 x i = T"^2 i=2 z% / ^- = ^(expl- 2 ) — z ~ 1); Si=2^A' 
X^f=i iz l /i\—z = zexp( 
and Eq. (JTSJ), we have: 



Si=i iz l /i\—z = zexp(z)z — z, Y^=2 x i hi^i^/T) = T ln 2; ^=2 = ^ lnz(exp(,z) — 1), 



2 — ^72 

zTe z — Tz + ztanh(7T)/7 = — - — , 

1-m (l-m)7 1-m 
<t>(T, m, 7) = — ^— ln — 

+ lncosh(7T) + ztanh(7T)[l — ln^/7 

+T[e z (l- zlnz) - z(l-\nz) -1]. (21) 



Let us now consider an ansatz for < Si\e Ht \Si >: 

< Si\ exp[AN{T - T )]\Si >< S 1 \e~ HM ff t °\S l > 

= exp{N[A(T - T ) + 0(T o ,m, 7 )]}. (22) 

While calculating this expression via saddle point, we first find the extremal point T = 
e~ 7 to/N from the saddle point condition: 

A-*®. (23) 

The transition period t\ = Ne J Ti is defined from the condition, that the contribution < 
Si\e~ Ht \Si > into Z of Eq. (7) is larger than the contributions of other terms < Sj\e~ Ht \Si > 
(equal to e* according to Eq. (13)): 

exp(iV[0(To,m,7) + A{T X - T )}) > expiNe^T,), 
m A m (b(Tn,m, 7) 
A — e 1 A — e 1 

Thus Eqs. (21), (23)-(24) give the relaxation period T\ = e'^ti/N under the constraint of 
Eq. (JUD for the fitness A. 

There are several phases in dynamics. For < t < to, there is a random drift to the peak 
configuration Si. For t < t < ti, there is a growth in the value of pi, but the macroscopic 
majority is still out of the peak configuration. For t > ti, the macroscopic majority is near 
the peak configuration. 

Let us give an explicit expressions for the case 

7(1 — m) 



A 



< 1. (25) 



This is a typical biological situation for observing l-m<l. In this case, as we can check 
later, J 1 ~ (1 - m) < 1, thus one can replace ztanh(7T)/7 — ► zT and derive a simplified 
system of equations: 

1 YY) 1 YJ1 

0(T,m, 7 ) = ^-[ln 7 ^--(l + lnz)]+T(e 2 -l), 

m , 1 — m 

Tze z = 

(26) 
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Then To = (1 — m)/[2(l + A) ln(l + A)]. Thus for the relaxation period t = Tie 7 iV, one has 
an expression: 

ln 2eln(A+l) 

t x = (l-m)N —S^IUh (27) 
1 v ; 2(Ae-"i - 1) y J 

Equation (|2*7J) gives relaxation period from the original distribution, concentrated at the 

configuration with the overlap Nm with the peak fitness configuration, and mutation per 

site 1 — q = i/N. The physical meaning of the term ^~^ N [ s trivial (for the case of 

infinite population): the relaxation period is proportional to the Hamming distance. We 

can understand also the term (Ae~ 7 — 1) in the dominator: it is a natural consequence of the 

fact that relaxation period should diverge at the error threshold Ae" 1 — > 1. Our derivation 

is valid when the condition of Eq. ()25|) is satisfied. An estimate for the t\ has been given in 

Refs. |l3, 



14. 

h = Ae- N ^-") - 1 = Ae-J-Y (28) 
We note that Eq. (|2*Hj) is qualitatively correct and consistent with Eq. (}2Tj) for the case 
N(l — m)/2 = 1 considered in that works. Our derivation is rigorous only for a large 
number of flipped spins, i.e. N(l — m)/2 >> 1. For a small number of flipped spins 
considered in Refs. jl3L Il4j , we still can not derive an exact analytical formula. 

Let us briefly consider a simple case of a dynamic fitness landscape: a fitness peak A(t) 
in the first configuration Si, which changes with the time. Now for the < Si\e~ Ht \Si >, we 
have exp[e~ 7 J * A{t)(Lt]. Equations ()23jl and ()24j) transform into 

A(t ) = d ^-, 0(T o ,m, 7 ) + [ Tl A(r)dr] > e 7 T\. (29) 

Now could be a very rich phase structure with different solutions for T . For the T\ = 
tie -7 /N, we have: 

A = £*gL, T 1= *T,-*e«™l. (30) 
Ti - T A - e 7 A - e 7 

Now A is replaced with a mean value. For the case of A ^> 7(1 — m), we again have Eq. 
(JZZj), only with A -> A. 

For A ^> 1, we can calculate the relaxation time from an original uniform distribu- 
tion on a static landscape: pi = 1/2 N . For this purpose, we compare the contribution 
< Si\e~ Ht \Si >= 2~ N expL4e~ 7 t] with exp(i) (sum of other contributions) for their contri- 



S 



butions to Z of Eq. (J7J) and find that 

AHn2 



Ae~"t - 1 



(31) 



To derive the steady state distributions of pi, we can set dpt/dt = in Eq. (JTJ). For 
A > 1 we can derive that pi = q N [(l — q)/q] d ^ 1 '^ and and the result obtained in Ref. j^: 

Let us briefly consider the case of two isolated flat peaks in fitness landscape with fitness 
heights Ai and A2, and widths g± and g%- The peak of height Ai has gi one-flip neighbors 
of the same height. A simple consideration gives for the effective fitness AAA + ^(1 — q)]. 
Thus the Svetina-Scuster phenomenon [r| for two peaks appears at Ai[l + gi(l — q)] = 

A 2 [l + (l-q)g 2 ]. 

In 1971, Eigen [lj found an exact error threshold for his model from information theory 
arguments. After more than 30 years of different approximate or numerical investigations 
of the Eigen model, we have found the exact dynamics of the model presented in Eqs. 
(21), (23), and (24). Our Eq. (27) gives the relaxation periods with a hi gh d egree of 



accuracy 0(1 — m) 2 ~ (d/N)2, it is more accurate than Eq.(28) derived in [13J, |l4j. In 
[l| we compared the accurate result of this work Eq.(27) with the corresponding relaxation 
period of parallel scheme to conclude, that even at the limit of vanishing mutation rates 
two mutation schemes give a finite (nonvansihing) difference in relaxation periods. Therefor 
there is at least one situation, when our exact Eq.(24) or accurate approximation Eq.(27) 
give new qualitative result. We have also applied the similar method to study a simple case 
of dynamical environments and obtained Eqs. (}2"9"j) and (j3*Uj) . The more involved situations 
with a very rich and interesting phase structure |l2| as well as the virus-immune system 
coevolution [14| can also be investigated by our method. The main open problem is an 
application of the same method to the finite population case. In this case the search of a 
a peak configuration could be exponentially large function of N, instead of a linear in Eq. 
(27). We hope that progress in this direction is possible in the near future, considering funnel 
like fitness landscapes. In any case in this work we considered Eigen model's dynamics as a 
some statistical mechanics problem and exactly solve it. 
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